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Abstract 

New model equations are derived for dynamics of self-aggregation of 
finite-size particles. Differences from standard Debye-Hiickel and 
Keller-Segel |2] models are: a) the mobility fx of particles depends on 
the locally-averaged particle density and b) linear diffusion acts on 
that locally-averaged particle density. The cases both with and with- 
out diffusion are considered here. Surprisingly, these simple modifica- 
tions of standard models allow progress in the analytical description of 
evolution as well as the complete analysis of stationary states. When 
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ji remains positive, the evolution of collapsed states in our model re- 
duces exactly to finite-dimensional dynamics of interacting particle 
clumps. Simulations show these collapsed (clumped) states emerg- 
ing from smooth initial conditions, even in one spatial dimension. If 
jj, vanishes for some averaged density, the evolution leads to sponta- 
neous formation of jammed patches (weak solution with density having 
compact support). Simulations confirm that a combination of these 
patches forms the final state for the system. 

Keywords: gradient flows, blow-up, chemotaxis, parabolic-elliptic 
system, singular solutions 
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1 Historical perspective of continuum models 
of self-aggregation 

Many fields of physics, chemistry and biology deal with the problem of de- 
scribing the evolution of the macroscopic density of a large number of parti- 
cles, which self-consistently attract each other over distances large compared 
to their mean separation. A related paradigm arises in biosciences, partic- 
ularly in chemotaxis: the study of the influence of chemical substances in 
the environment on the motion of mobile species which secrete these sub- 
stances. One of the most famous among such models is the Keller-Segel 
elliptic-parabolic system of partial differential equations j2] , which was intro- 
duced to explore the effects of nonlinear cross diffusion in the formation of 
aggregates and patterns by chemotaxis in the aggregation of the slime mold 
Dictyostelium discoidium. The Keller-Segel model consists of two strongly 
coupled reaction-diffusion equations, 

p t + divJ = J = p,u($)V$ - DVp e$ t + L$ = 7p , (1) 

expressing the coupled evolution of concentration of organisms (density) p 
and concentration of chemotactic agent (potential) $. These challenging 
nonlinear equations were posed for chemotaxis as an initial value problem 
with Neumann boundary conditions initial data. The constants 7, D, e > 
are assumed to be positive and J is the flux of organisms (particles). The 
linear operator L is taken to be positive and symmetric. For example, one 
may choose it to be the Laplacian L = —A or the Helmholtz operator, 
L — 1 — a 2 A. For the moment, the functional dependence of the mobility 
/i will be left unspecified. 1 The "simplified" Keller-Segel model is obtained 
by setting e = in (1). Because of its fascinating mathematical properties 
- such as finite-time concentration of the density p into Dirac delta-function 
singularities, starting from smooth initial conditions - and its relevance to 
chemotaxis, the Keller-Segel system has attracted sustained interest in math- 
ematical biology. To avoid excessive citation here, we shall only refer to the 
review article [3] where over 150 references on the role of in the chemotaxis 
problem may be found. 

1 One recovers the precise expression of Keller and Segel [5] by rewriting the mobility 
term as //(<&) V$ = Vx( ( E ) ), where x(^) is the "sensitivity function." The development 
below introduces a different functional dependence in the mobility fi. 
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The applications of equations in the form of system JI} in physics and 
chemistry may be traced back into the 19th century. The history of these 
equations is a story of recurrent rediscovery, each time as a leading order 
description of singularity formation, occurring as the tendency to diffuse is 
eventually overwhelmed by long-range forces of attraction. Quite likely, their 
story has not finished yet and these equations will be needed again for the 
insights they provide into nonlinear multiscale dynamics in the competition 
between coalescence and diffusion. 

Three main modeling steps figure in the derivation of and similar 
equations. First, the flux J must be modeled in terms of the local particle 
density, its spatial gradient and the gradient of the potential. Second, an 
equation must be formulated for determining the potential $ from the local 
density p. Third, conservation of mass, charge, or particle number must be 
introduced as the evolution equation for density. 

Historically, it seems that Debye and Hiickel in 1923 were the first to put 
all three of these modeling steps together. They derived the evolutionary 
system ([TJ in their article pQ on the theory of electrolytes. In particular, the 
simplified model with e = in may be found as equations (2) and (2 ') in 
Debye and Hiickel £Q. Consequently, the simplified evolutionary system (0) 
with e = may also be called the Debye-Hiickel equations. 

The Debye-Hiickel equations (JTJ) with e = also appeared very early 
in astrophysics as the Poisson-Smoluchowski system, a particular Fokker- 
Planck equation describing a self-attracting reinforced random walk of parti- 
cles (with mass, but without inertia) diffusing by local collisions and drifting 
against friction under their mutual long-range attraction. As is well known, 
gravitational force leads to collapse in this case, when mutual attraction fi- 
nally prevails over diffusion and friction. See Chandrasekhar (1939) [3] for 
an exhaustive historical review from the viewpoint of stellar formation. 

Interest in the Debye-Hiickel system was recently revived, when the "Nernst- 
Planck" (NP) equations in the same form as (0) re-emerged in the biophysics 
community, for example, in the study of ion transport in biological channels. 
In this case, the flux J is called the Nernst-Planck (NP) particle flux for the 
ion current density J depending linearly on the gradients (Vp, V$) and L 
is the Laplace operator. See Barcilon et al. (1992) j3] and Syganow and 
von Kitzing (1999) [7j for recent reviews of ion transport in biological media. 
The same elliptic-parabolic system had also surfaced earlier as the "drift- 
diffusion" equations in the semiconductor device design literature, Selberherr 
(1984) jS]. Similar elliptic-parabolic equations (with hyper-diffusion, instead 
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of diffusion) also arose recently as a leading order description of molecular 
beam epitaxy j^j. A variant of the system re-appeared even more recently 
as a model of self-assembly of particles at nano-scales [TTj . 

2 Derivation of aggregation equation for 
finite-size particles and variable mobility 

In this paper, we modify the class of Debye-Hiickel equations (0) with e = 0, 
as a model of the aggregation of interacting particles of finite size. This prob- 
lem is motivated by recent experiments using self-assembly of nano-particles 
in the construction of nano-scale devices J2|-[2n|- A colloidal solution of 
50nm-size particles is deposited on a grooved substrate; evaporation and re- 
ceding contact line drags the particles into the channels. It has been shown 
[TT] that the most important phenomenon effecting the distribution of parti- 
cles in the channels is the capillary interaction between the particles at last 
stage of water evaporation, when all water and particles are confined to the 
channels. Examples of particle self-assembly in nano-channels is shown inFig. 
1. Fundamental principles underlying the mathematics of self-assembly at 
the nano-scales are non-local particle interaction and nonlinear motion due 
to variations of mobility at these scales. 

Let us start by deriving the effective interaction potential between nano- 
particles. This will reveal the length scales involved in the problem and 
characterize the type of interaction between the particles. Due to the effects 
of surface tension, the air/water interface is deformed near a particle. On 
the other hand, water is attracted to the substrate surface by van der Waals 
forces. Denote the thickness of the undisturbed water layer by Hq and let h 
be the elevation from this equilibrium level. If we call the potential of van der 
Waals interaction U v dw{H), then for small deviations from the equilibrium 
surface the balance of surface tension and van der Waals force gives 

7A/1 = F h, F = -^^(H = H ) (2) 
A typical expression for van der Waals potential includes three terms [21] 

UvMH) = ^-§- K -Ce~^ (3) 

with constants A, B, C, with k > 3 and Iq being the interaction length. A 
numerical estimate for these values from [21] shows that for H ~ 50 nm, 
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the first term on the right-hand side of is larger than the second and 
third terms by factors 10 9 and 10 18 , respectively. Here, A is called Hamaker 
constant, which assumes a certain value for a given pair of materials. 2 Then, 
(j2J) can be re-written in terms of capillary length l c as 



where we have applied the term capillary length even though the stabilizing 
potential is due to van der Waals forces and not gravity. Numerical values for 
Hamaker constant for water on silicon is about 30zJ = 30 • 10~ 21 J, surface 
tension of water/air interface is 0.07N/m. Consequently, for Hq = 50nm, 
equation (j3J) yields capillary length l c ~ 300nm, or several times the particle 
diameter. 

It is now clear that in spite of very different physics and length scales, 
the phenomenon of self-attraction of particles at nano-scales bears not just 
qualitantive, but also quantitative resemblance to the well-known Cheerios 
effect: self-attraction of floating bodies due to surface tension. For floating 
cheerios, the typical length of interaction is the usual capillary length (several 
mm), which is also of the order of particle size. It has been shown [221122] that 
two particles separated by a distance I have interaction potential proportional 
to exp(— l/l c ) in one dimension, and K (l/l c ) in two dimensions (Ko(x) being 
the modified Bessel function of the first kind). The force between the particles 
is then proportional to Ki(l/l c ), and this result holds for both floating and 
partially submerged particles. Mathematically, these interaction potentials 
arise from inversion of Helmholtz operator in (@J). 

These physical considerations show that interaction potentials propor- 
tional to the Green's function for the Helmholtz operator plays a fundamental 
role in particle self-assembly across surprisingly many orders of magnitude. 
Thus, we shall consistently use this interaction potential in all our numerical 
simulations. Nevertheless, we shall keep our framework general enough that 
our method of modeling nonlocal interactions among different particles may 
be extended to the other physical problems of interest described in the pre- 
vious section. Our method takes into account the change of mobility due to 
the finite size of particles and the nonlocal interaction among the particles. 
The local density (concentration) of particles is denoted by p. Suppose the 

2 Usually, the numerical vallue of A is given in zJoules, with z being used for zepto= 



l 2 c Ah -h = 
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particles interact pairwise via the potential — G(|r|). The total potential at 
a point r is 



where * denotes convolution. (The minus sign is chosen so G > for attract- 
ing particles.) The velocity of the particle is assumed to be proportional to 
the gradient of the potential times the mobility of a particle, p. The mobility 
can be computed explicitly for a single particle moving in an infinite fluid. 
However, when several particles are present, especially in highly dense states, 
the mobility may be hampered by the particle interactions. These consider- 
ations are confirmed, for example, by the observation that the viscosity of 
a dense suspension of hard spheres in water diverges, when the density of 
spheres tends to its maximum value. Many authors have tried to incorpo- 
rate the dependence of mobility on local density by putting p = p(p) 
It is common to assume that the mobility is a function of density, which 
tends to zero when the density tends to some maximum value, assumed to 
be Pmax = 1, *-e., jti(pmax) = 0. Vanishing mobility leads to the appearance 
of weak solutions in the equations, to singularities and, in general, to massive 
complications and difficulties in both theoretical analysis and computational 
simulations of the equations. In contrast, the Keller-Segel model specifies 
so the mobility is taken in that function of the concentration 

of chemotactic agent (potential) $. 

The goal of this paper is to suggest and investigate the implications of 
the idea that the mobility p should depend on the averaged density p, rather 
than either the potential, or the exact value of the density at a point. This 
assumption makes sense from the viewpoints of both physics and mathemat- 
ics. From the physical point of view, the mobility of a finite-size particle must 
depend on the configuration of particles in its vicinity. While attempts have 
been made to approximate this dependence by using derivatives of the local 
density, it is much more natural, in our opinion, to assume that the local 
mobility depends on an integral quantity p which is computed from the den- 
sity as p = H * p. Here H(y) is some filter function having, in general, short 
range compared to the potential G. Several filter functions are possible, with 
examples being H (r) = 5(r) (5-function), H(r) = exp(— |r|//)/(2Z) (exponen- 
tial, or inverse-Helmholtz in ID) or, in d dimensions, H(r) = 0(\r\ —l)/(2l) d . 
The last is the top-hat function, whose value is unity when |r| is between — / 
and I, and vanishes, elsewhere. The normalizing factor of (2l) d (d being the 
dimension of space) is introduced so that f_™ H(r)dr = 1. Alternatively, 
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one may assume a filter function H(r) with zero average: f_™ H(r)dr = 0. 
The latter assumption may be useful in crystal growth models, where the 
mobility does not depend on the absolute level of the material. Rather, in 
such models the mobility depends on the relative positioning of particles. 
The mathematical analysis and the reduction property derived in this paper 
hold, regardless of the shape of the filter function H and the potential G is, 
so long as they are both nice (e.g., piecewise smooth) functions. 

Aim of the paper The object of this paper is the analysis of the following 
continuity equation for density evolution, 

dp 

-rf = -divJ with J = - D Vp - p(p)pV$ . (6) 
dt / s v ' 

Continuity equation Particle flux 

Here D is a constant diffusivity, while p and <3> are defined by two convolutions 
involving, respectively, the filter function H and the potential G, 

p = H * p and $ = G * p . (7) 

While this equation preserves the sign of p, we shall see that it allows the 
formation of 5-function singularities even when D > and in the case of one 
spatial dimension. 

Alternatively, system (|6I7|) with D = can be represented geometrically 

as 

^(pdVol) = along ^ = J(p) , 

where the vector field corresponding to the current J(p) is defined in terms 
of the pointwise density p > and two prescribed functionals of p: an energy 
E(p) and a mobility p(p) > 0, as 

3(p) = -pp(p)V—. 

dp 

Here p = /J*p>0isan average density. As a result, 

d r 1 



4-E(p) = - ( —- |J(p)| 2 dVol < 
dt y 1 J pu(p) 1 v n 



Thus, the flow dx./dt = J causes the energy functional E(p) to decrease to- 
ward its minimum value and the flow is Lyapunov stable, provided the energy 
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functional E(p) is sign definite and p(p) is isolated from 0. We may regard 
J as a vector field defining an infinitesimal action of the diffeomorphisms 
on the space of positive maps p acting on a manifold M.. Remarkably, this 
action produces singularities in which the pointwise density concentrates in 
finite time on subspaces embedded in the manifold M.. 

Subcases The generality and predictive power of the model ()6I7|) can be 

demonstrated by enumerating a few of its subcases, as follows. When H(r) = 
5(r) and G(x) = e~' r '^ the system reduces to the generalized chemotaxis 
equation [21]. For the choice G = S'(r), H = 5'(r) one obtains a modification 
of the inviscid Villain model for MBE evolution [§| (with extra factor of p 
in the flux). If the range of G (denoted by I) is sufficiently small, we can 
approximate (at least formally) the integral operator in $ = G * p as a 
differential operator acting on the density p, namely, 

$ = p + /V • (IV p) . 

As was demonstrated in [11] . equation (JBJ) then becomes a generalization 
of the viscous Cahn-Hilliard equation, describing aggregation of domains 
of different alloys. In the case that p, is a constant, H(x — y) — 8(x — 
y) , and G is the Poisson kernel, then equation © becomes the drift limit 
of the Poisson-Smoluchowski equation for the interaction of gravitationally 
attracting particles under Brownian motion jl] . 

All the particular cases described above require a singular choice of the 
functions G and H. However, the generalized functions required in these 
cases may be approximated with arbitrary accuracy by using sequences of 
nice (for example, piecewise smooth) functions. We shall concentrate on 
cases where the functions H and G remain nice, and derive the results in 
this more general, more regular, setting. Thus, regularization in this endeavor 
introduces additional generality, which enables analytical progress and makes 
numerical solution of the equations easier. 
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3 Evolution of density as mass conserving gra- 
dient flow: the nonlocal Darcy's law and 
energetic considerations 

Let us have another look at our motivation of equations (J6I7|) . this time mak- 
ing a connection with gradient flows and thermodynamics. We will show that 
there is a naturally defined free energy which remains finite, even when the 
solutions cease to be smooth and are replaced by a set of delta-functions. 
We show that the energy remains finite for all choices of p(j>). In contrast, 
the traditional approach considering the dependence of mobility on the un- 
smoothed density p(p) would fail by allowing the energy to become infinite. 
This additional energetic argument reinforces the choice of the regularized 
model in (J6I7|) . The calculation will be performed in arbitrary spatial dimen- 
sions. 

For the energy to be well-defined, we require that the function G describ- 
ing interaction between two particles is everywhere positive, so the interac- 
tion is always attractive. In one dimension, we also require that the kernel 
function G is symmetric and in n dimensions that the interaction is central. 
This assumption is physically viable for all the historical cases described in 
the introduction. Our derivation for the free energy will remain valid for an 
arbitrary filter function H . 

We start with equation © for mass conservation, in the case when the 
mobility p(~p) in the particle flux J is not constant. We shall derive a variant 
of equation © as a gradient flow in the sense that, 

| = -grad£| p , or (| , = - , 0} , 

where (/ , <ft) = J fcf)d n x is the L 2 pairing, for a suitable test function <fi and 
where E is the following energy, 

E[p] = L>/p(logp- l)dx + ~Jp<&dx with $ = G*p. (8) 

The first term in this expression for the energy E[p] decreases monotonically 
in time for linear diffusion. The second term defines the if -1 norm for 
the choice G(x) = exp(— It may be regarded as a generalization of 
this norm for an arbitrary (but positive and symmetric) function G. Such 
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negative Sobolev norms remain finite, even when the solution for the density 
p concentrates into a set of delta functions. 

The variation 5E/5p of the free energy E[p] in equation (JSJ) is 

5E[p}= J(D\ogp + $)5pd n x 1 

where we understand dp as arising from the flow of a diffeomorphism, as for 
example in Namely, we specify 

6p = -div(p/i(p)w) , 

for an arbitrary variation of p determined by the function which is assumed 
to be smooth. 3 The corresponding variational derivative is given by, cf. 

SE[p] = (-^, tf) = - J(Dlogp + $)div(pp(p)W)d n x 

= - J ^dw(p(p){DVp + pV<S>))d n x. (9) 

Consequently, the free energy E[p] in equation (jHJ) produces the following 
gradient flow: 




) = (div(/i(p)(DVp + P V$)),*). 



The resulting variant of equation (jUJ) is 

^ = -divJ with J = - p(p)(DVp + pV®)) , (10) 

vi2 , ' V v ' 

Continuity eqn Modified particle flux 

in which the linear diffusivity of density p is modified by the nonlocal mobility 

3 This specification of the density variation formally requires the solution to remain 
diffcrentiable and the product pp,(fi) not to vanish. In principle, these conditions would 
be violated by the formation of weak solutions, in which the density concentrates into 
delta- functions. However, we will verify a posteriori that the gradient flow properties of 
the smooth solutions discussed in this section are also preserved for the weak solutions. 
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Energetics The free energy E[p] in JHJ) decreases monotonically in time 
under the evolution equation (fTU|). A direct calculation yields, 



- [(Dlogp + <5>)dwJd n x = - [ ——r \3\ 2 d n x. (11) 

J J DUiP) 



dE[p] 

dt J y '~ or '->■■- J p ^ ( p) 

According to this equation, provided pp(p) > 0, the rate of decay of free 
energy E[p] given in (jSJ) defines a Riemannian metric in the particle flux, 
cf. J2S1- We shall see that when D = 0, the resulting conservative motion 
of finite-size particles drifting along the gradient of $ leads to a type of 
'clumping' of the density p into a set of delta functions. One may check 
that this monotonic decrease of energy persists for more general functions, 
including weak solutions supported on 5-functions, as discussed below. 



4 Non-vanishing mobility: weak solutions 
(clumpons) and their role in long-term dy- 
namics in ID 

Everywhere in this paper, we shall assume that the potential is purely at- 
tractive, so G(x) > for all x. This assumption is used for two main reasons. 
First, it suits the physics of the motivating problem, namely, mutual attrac- 
tion of nano-particles ^Uj- Second, restricting to purely attractive interac- 
tions allows one to separate all equations of the type ()6I7|) into two different 
classes. The physical property separating the two classes is whether the mo- 
bility p(p) is strictly isolated from zero, i.e., p(p) > p > or one may allow 
p(p ) = for some p = p . In the first case, there is nothing to resist the 
mutual attraction of the particles and the final state of the system is a single 
5-function no matter what the form of the mobility dependence, as long as 
mobility is strictly isolated from zero. In the second case, vanishing of the 
mobility at some maximum density blocks the motion when the densities be- 
come too large and prevents total collapse. Instead of collapse, this produces 
isolated patches of solutions of constant density. It is interesting that also in 
this second case, the stationary solutions remain exactly the same regardless 
of the dependence of p on p, and we can describe these stationary solutions 
analytically. Moreover, we shall demonstrate that these solutions are stable 
and any initial condition separates into a set of isolated stationary patches 
of this type. 
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This section is devoted to analytical description of dynamics in the case 
when /z(p) being strictly isolated from zero. In simulations, we have assumed 
p(p) = 1 for simplicity, but all our results will remain valid for arbitrary 
dependence of /i on p, as long as p always remains positive. 



4.1 Formal weak solution anzatz 

This section considers motion under equations (|6I7|) in one spatial dimension 
for the case of vanishing linear diffusivity, D = 0. Substituting the following 
singular solution ansatz 

N N 

p(x,t) = $>i(*)*(z -*(*)) , p{x,t) = y £w j (t)H(x-q j {t)), (12) 
i=i j=i 

into the one- dimensional version of equation JBJ) with D = and integrating 
the result against a smooth test function yields 



J <j> Pt~ (pp{p){G * p)^j x dx = J (fi(x)^Twi5(x - qi)dx 

i=i 

. N N 

+ / $\ x )^ w i(qi + ^ w j(t)fi(p)G'(x - qj))5(x - qi) dx 



Hence, one obtains a closed set of equations for the parameters Wi(t) and 
qi(t), i = 1, 2, . . . , N, of the solution ansatz (|12|). in the form 

N 

Wi(t) = 0, qiit) = -J2 w jP(p) G '(qi- qj) (13) 

where 

N 

P = w m H(q m (t)) (14) 

771=1 

Thus, the density weights Wi(t) = Wi(0) = Wi are preserved, and the posi- 
tions qi(t) follow the characteristics of the velocity u = —p(j))VG * p along 
the Lagrangian trajectories at x — qi(t). This result holds in any number 
of dimensions, modulo changes to allow singular solutions supported along 
moving curves in 2D and moving surfaces in 3D. Fig. 2, demonstrates that 
the solutions (|13|) do indeed appear spontaneously in a numerical simulation 
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of equation (JBJ) with D = 0.02. Hence, they are ubiquitous and dominate its 
dynamics. The simulation started with a smooth (Gaussian) initial condi- 
tion for density p. Almost immediately, one observes the formation of several 
singular clumpons, which evolve to collapse eventually into a single clumpon. 
Observe that the mass of each individual clumpon remains almost exactly 
constant in the simulations, as required by equation (|13j) . Also note that the 
masses of two individual clumpons add when they collide and "clump" to- 
gether. Eventually, all the mass becomes concentrated into a single clumpon, 
whose mass (amplitude) is exactly the total mass of the initial condition. 

4.2 Energy decay close to the final state and estimates 
for collapse time 

Normally, systems approaching an equilibrium state tend to evolve more 
slowly as they approach the equilibrium. If the rate of approach diminishes 
linearly with the distance from equilibrium (in some sense), for example, 
one obtains an exponential decay towards the equilibrium. Alternatively, for 
finite-time singularities, the rate tends to diverge to infinity in a power-law 
fashion. In contrast to these familiar examples, the system (1617)) approaches 
the singularity at a constant rate. That is, the rate of approach to singularity 
never diverges. Consequently, one may predict the formation of singularities 
and even predict their evolution after they have formed. 

This surprising result may be demonstrated by deriving an alternative 
form of energy dissipation. Direct substitution of a single 5-function for 
density into (fTTj) is mathematically ambiguous, as it would lead to improper 
operations with ^-functions. Instead, let us notice that the evolution of the 
energy E describing the gradient flow of (Jb17j) can be expressed as follows 

— — = — / G(x — x )p(x — - — dxdx = 
alt J y IHy ' dt 

+ / $div (p/i(p)grad$) = - / p/i(p) ( V$) 2 (15) 

Note the difference between this formula and (fTTj) . Formally, (|TK|) is the same 
as 1)1 but we can substitute the delta- function ansatz for weak solutions 
directly into (JT3J). 

Numerical simulations, energetic considerations and physical intuition 
suggest that the final state of the system with strictly positive mobility 
Mp) — > should be a single clumpon. One may ask how the time 
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necessary to collapse to this final state varies with the initial conditions. 
An exact analytical answer to this question seems out of reach and we shall 
provide a numerical simulation for various initial conditions. However, a sur- 
prisingly simple and accurate analytical estimate can be made here. To start, 
let us compute the rate of energy dissipated by a single clumpon, which is the 
final state of our system. We assume p = M5(x), where M is the total mass. 
Assuming that H(x) is regular, H * p(x) = MH(x). Let us also assume for 
simplicity that G{x) = Goexp(— |x|/a). Direct computation gives 



dE M 3 
— \p = clumpon] = —Mfi (MH(0)) (V$(0)) 2 = -p (MH(0)) Gl (16) 

dt or ' 



We may now estimate how long it will take for an arbitrary system to collapse 
into a single clumpon. If the initial energy of the system is E , and initial 
mass is M, the energy of a single clumpon will be 

E f = M 2 G . (17) 

Then, the approximate time to collapse to a single clumpon from any initial 
conditions is given by combining and (|17|) : 

E - E f E - M 2 G 

* dE/dt[p = clumpon] M 3 /i (MH(0)) G 2 [ ' 

In derivation of (|18|) we have assumed that dE/dt does not change much along 
the trajectories, so the initial conditions are close enough to a clumpon. Thus, 
this estimate for the time of collapse depends on only two integral quantities: 
initial energy and mass for initial conditions close to the final state. 



4.3 Rate of blow up 

The following analysis of the evolution of a density maximum reveals that the 
clumping process results from the nonlinear instability of the gradient flow in 
equation (jUJ) when D = 0. For a particular case G = G e~ N/Q and p(p) = 1, 
one may show that for a high enough peak, a density maximum p m (t) = 
p(x m (t),t) becomes infinite in finite time. The motion of the maximum is 
governed by 

d 1 1 

-rPm = — (p 2 m - PmHXmj) > — (fa - p m M) , (19) 



15 



where M = J pdx is total mass and we have used the fact that $ satisfies 
$ — a 2 § xx = p. The last inequality holds, because G < 1 is bounded and p 
is everywhere positive. Thus, if at any point the maximum of p exceeds the 
(scaled) value of the total mass, then the value of the density maximum p m (t) 
must diverge in finite time. This divergence produces (^-functions in finite 
time. From (|19j). the density amplitude must diverge as p m ~ a 2 /(t — t). 
The formation of singularities in Fig. 2 occurs both at the maximum, and 
elsewhere. The subsidiary peaks eventually collapse with the main peak. 

4.4 Dynamics of inflection points and collapse 

Let us discuss in more detail the process by which density singularities are 
formed. Suppose the the initial condition contains only one density maxi- 
mum. Then, one would expect the first singularity to form at the position 
of this maximum. As this density singularity forms, mass flows toward it. 
This causes a local reduction of density in the neighboring region from which 
the mass is flowing. Whenever this local reduction of density due to the 
formation of the singularity is sufficient to form a new maximum in density 
away from it, then the process of singularity formation starts again there, 
and so forth. In principle, this process could form an infinite number of 
density singularities. However, we propose two reasons why the expected 
number should be finite. First, the range of interaction is set by the length 
scale in H. Thus, one might expect no more clumpons to form than the 
ratio of the interaction range to the domain size. Second, in the formation 
of singularities from two nearby maxima, one of them may become stronger 
than the other and entrain it, thereby suppressing the formation of the sec- 
ond singularity. In practice, one sees the formation of considerably fewer 
clumpons than the number estimated by the ratio of the interaction range 
to the domain size. Hence, one may conjecture that the formation of singu- 
larities does involve some competition between neighboring density maxima. 
This conjecture could be tested by studying evolution from initial conditions 
containing several density maxima separated by varying distances which are 
comparable to, or smaller, than the average length scale in H. 

A quantitative evaluation of this heuristic argument may be obtained 
by computing the position of the inflection point for the averaged density p. 
Although no analytical formula for the motion of inflection point is available, 
one sees numerically that the inflection point quickly converges to the position 
of the singularity. The evolution of the inflection point corresponding to the 



16 



solution from Fig. 2 is shown in Fig. 3. 

5 Formation of jammed states when mobility 
fi approaches zero 

5.1 Competing length scales of H 8z G 

An interesting limiting case arises when the scale of non-locality of H is much 
shorter than the range of the potential G. Formally, this limit corresponds 
to H(x) — > 5(x). In practice, we may select a sequence of piecewise smooth 
functions H e (x) = exp(— |x|/e)/(2e) which converge weakly to a 5-function. 
For each function H e (x) in this sequence, no matter how small (but positive) 
the value of e, the exact ODE reduction (|13j) still holds. So, what happens 
in the limit of very small epsilon, as e — > 0? In investigating this limit, 
we performed a sequence of numerical simulations for a fixed choice of the 
function G{x) = exp(— \x\/a) while varying H e (x) = exp(— x/e)/(2e) over 
a sequence of decreasing values of the ratio e/a. This simulation is shown 
in Fig. 2 for e/a = 1/10 and it demonstrates the formation of flat clumps 
of solutions. The mechanism for this phenomenon is the following. The 
vanishing mobility at p — 1 caps the maximal density at p — 1 in the long- 
term. This leads to the appearance of flat mesas in p(x) for large t. On the 
other hand, when H(x) — > S(x), one finds ~p(x,t) — > p(x,t) pointwise in x, 
which forces p(x, t) to develop a flat mesa, or plateau, structure in which the 
maximum is very close to unity, as well. This is precisely what is predicted 
for the chemotaxis equation with H(x) = S(x) and p(p) = 1 — p |24j . 

In the limit H —> 5, model ()6I7|) recovers ordinary diffusion of local den- 
sity. This is a singular limit, because it increases the order of the differenti- 
ation in the equation. Since ordinary diffusion is known to prevent collapse 
in one dimension [Hj, this singular limit should be of considerable interest for 
further analysis. 

5.2 Stationary states 

Numerical simulations of time-dependent problem (|6I7|) show that the solu- 
tion converges to well-defined states with flat mesa peaks when n(p) reaches 
for some value of p. In the remainder of this section, we shall concentrate 
on the analytical description of the evolution in this case. For simplicity of 
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formulas, we shall assume that the critical value of p is normalized to be 1, 
i.e., /i(l) = 0. In simulations, we shall take p = 1 — p. However, all theoret- 
ical results, in particular, exact expressions for stationary sates, remain true 
for arbitrary dependence of p(p), as long as p(p) reaches at some value of 

We distinguish two classes of stationary states, which differ mathemati- 
cally and physically. Both classes of stationary states will describe a clump 
of particles whose density p(x) is a generalized function with compact sup- 
port, which we will call a patch. The difference between these classes lies 
in the physical reason for the vanishing of local velocity. Remember that 
local velocity is written as u = pfflS/Q. The solution is stationary if u 
vanishes at every point inside the patch. This can be achieved in two ways. 
First, V$ may be exactly zero everywhere inside the patch. In this case, 
the net force acting on every particle vanishes. These solutions are called 
equilibrium states. Second, it may also happen that the mobility p vanishes 
at every point inside the patch. Physically, this corresponds to a clump at 
maximum density whose motion is prohibited, although the net force on each 
particle may not be zero. These solutions are called jammed states. Let us 
now study both the equilibrium and jammed states in more details. As we 
shall see, an analytical solution describing each of these states can be found. 
These solutions will also elucidate the nature of competing length scales in 
G and H . The stability of these solutions and their selection mechanisms 
will also also be demonstrated. 

5.3 Equilibrium states 

Let us assume for now that both the potential attraction G(x) and the filter 
function H(x) are proportional to the Green's function for the Helmholtz 
operator. This is both the case we use in numerical simulations, and also, 
incidentally, the only case we have found that admits analytical solution for 
stationary states. Thus, we posit 

G(«)=«p(-M) tf W = _L eX p(-M) (20) 

We are looking for stationary states p{x) that are weak solutions with com- 
pact support, for which p{x) = if \x\ > L. While it is impossible to find 
a smooth (or even piecewise smooth) function p(x) which is a stationary 
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solution of (jHJ), one can find a stationary solution of the following form 

p{x) = w5(x + L) + w5(x - L) + ( °' } x j > L T (21) 



where the value of the constants w and L are to be determined. 

For (|21|) to be a stationary (weak) solution, the following conditions must 
be satisfied: 

= (G * p)(x) = const, \x\ < L , . 

p{x) = (H *p)(x) = 1, x = ±L. ^ 

Substitution of (|2T|) into the first condition of (j22l gives 

= wG(x -L)+ wG{x + L) + J G(x - y)p (y)dy = const (23) 

It is crucial to notice that if G(x) is given by (120(1 . then the integral in (|2*3*|) 
must yield a function proportional to a sum of exponentials, which is only 
possible if po(x) is a constant: po(x) = p . Performing the integral gives 



$( x ) = we ~\ x ~ L \' a + we~ lx+Ll/a + p « (2 - e^ x ~ Ll/a - e ^ x+Ll/a ) = const 

(24) 

which requires 

w = poa. (25) 

Enforcing the second condition in ()22|) will determine po as a function of L. 
As in (123(1 . we find an analytical expression for p 

= ^ ( e -l*-L|//J + e -|x + L|/^ + £ ( 2 - e-l-^ - e-l^ L l^) (26) 

Since ()26|) assumes identical values at x — ±L, we need only check that 
p(x — L) — 1. Substitution of x = L into (|2^j) gives 

Po [a (1 + exp (-2L//3)) + (3 (1 - exp (-2L//3))] = 2/3 (27) 

In principle, equation ()27() already determines the density inside the patch as 
a function of length. It is, however, more practical to determine the length 
as a function of total mass of the clumpon M. The mass contained in the 
solution (|2"T|) is 

M = p L + 2w = p (L + 2a). 
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Thus, we find an implicit condition for length L as a function of patch mass 
M: 

M[a(l + exp (-2L//3)) +(3(1- exp (-2L//3))] = 2/3(L + 2a) (28) 

These solutions are shown on Fig. 5. It is interesting to note that (J28)) de- 
fines a non-negative length L only if the mass M is sufficiently large. This 
is perfectly physical: if the mass is small, there is no reason for the solution 
to "spread out", so a single 5-function is produced. Incidentally, this will 
happen if the mass of the clumpon M = J pdx is such that max x p < 1, 
i.e., M < 2/3. Physically, we expect that equilibrium solutions are unstable. 
Indeed, imagine a set of particles on a line positioned at an equal distance 
from each other. If the force between each pair of particles is repulsive, this 
configuration is stable, and the real part of all eigenvalues of the linearized 
problem is negative. However, if the force between the particles changes to 
purely attractive without changing the particle positions, this leads to the 
sign change of the eigenvalues which correspond to a stable situation. Thus, 
these equilibrium states with fixed gaps between the particles "held together" 
by purely attractive forces must be unstable, and must exhibit large values 
for the real parts of the unstable eigenvalues. Numerical analysis of linear 
stability of equilibrium states in the continuum description confirms this in- 
tuitive physical picture. Namely, there exists a set of eigenvalues with large 
positive real parts. In addition, a fully nonlinear time-dependent simula- 
tion of ()6I7|) starting with initial conditions corresponding to even slightly 
perturbed equilibrium states shows very rapid deviation away from equilib- 
rium. Therefore, equilibrium states are unstable and will never be realized 
in nature. 



5.4 Jammed states 

The derivation of the jammed states is rather similar to that for the equi- 
librium states, so only a brief discussion of it will be provided. Let us again 
assume that stationary states p(x) are weak solutions with compact support, 
and p(x) — if \x\ > L: 



p(x)=wS(x + L)+w6(x-L) + { ' , (29) 



fo, 


\x\ 


1 Po(x), 


\x\ 



Now, for solution (}2*9*|) to be a jammed solution, we need just one condition: 
p( x ) = (H* p){x) = 1, |x| < ±L. (30) 
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Substituting (J23j) into condition (|3(J|) gives 

p(x) = wH(x -L) + wH(x + L)+ J H(x - y)p (y)dy = 1 (31) 

Again, it is essential that H(x) is given by (|2T)|) . so the integral in (J2Hj) must 
be proportional to a sum of exponentials, which is only possible if po(x) is a 
constant: po(x) = p . Equation (pTTj) transforms to 

= ™ L-\*-L\/P + e -|*+L|/^ + ^ ( 2 - e-l-^l// 3 - e-l^l/' 3 ) = 1 (32) 
which requires 

w = po = 1 (33) 

The mass contained in the solution (}2*T]) is 

M = p L + 2w = L + 2(3. 

The jammed states are illustrated on Fig. 6. 

To understand the nonlinear stability of jammed states, let us again ap- 
peal to the physical picture of particles on a line. Jammed states correspond 
to a set of finite-size particles pressed tightly together. Physical intuition tells 
us that such a state should be favorable for purely attractive forces between 
pairs of particles. Numerical analysis of linear stability of jammed states in 
the continuum approximation confirms this: the real part of the spectrum 
is isolated from by —D, where A = —D is the limiting point of spectral 
sequence. Thus, D should determine the time scale for convergence to sta- 
tionary solutions. That is, the rate of convergence to stationary solution is 
given by the time scale r ~ l/D. 

To verify these predictions, a fully nonlinear simulation starting with a 
smooth initial conditions in the form of Gaussian peak has been performed. 
The results of this simulation are given in Fig. 7. For D = 0.01 used in the 
simulation, the solution becomes practically stationary after t ~ 100. The 
small "bumps" at the position of the 5-functions in density correspond to a 
mismatch in the stationary solutions, which were derived for D = 0. The 
amplitudes of the "bumps" tend to zero when D becomes very small. 
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6 Jammed stationary states in two dimen- 
sions 



6.1 Exact jammed states in two dimensions 

While this paper focuses primarily on the evolution of particles in one di- 
mension, the real-world technological importance of self-assembly processes 
requires us to derive and analyze stationary solutions of our model in two 
dimensions. Numerical simulations of particle dynamics (21] show that there 
is a tendency for the formation of isolated clumps of either roughly circular 
shape, or, for high particle densities, roughly circular or elliptical voids in 
fully dense areas. To model these results, let us assume that the mobility 
p(p) vanishes at critical density p = 1, and p and p are connected through 
the two-dimensional Helmholtz operator 

p - p 2 Ap = p , so that p = H* p. (34) 

Inspired by these results we use the intuition developed in one dimension 
to seek jammed solutions in two dimensions which are fully dense p = 1 
inside an area D with additional distributed ^-function for density on the 
boundary dD, and p = in the exterior of D. Remarkably, when H is the 
Green's function for the Helmholtz operator, an exact analytical expression 
for several possible shapes of D can be found. Again, we shall emphasize 
that our results will hold for an arbitrary functional dependence p(p) as long 
as the mobility vanishes for some value of p = p*. For convenience, we have 
again rescaled the critical value to be p* = 1. 

Consider an orthogonal coordinate system (£, ri) for which the solution 
of Helmholtz equation separates variables. Suppose, in addition, that the 
boundary of D corresponds to one of the coordinate lines £ = £o ; and interior 
of D is given by £ < £o Then, we may seek the solution in the exterior of D 
in the form p(£, 77) = F(£//3)G(r]/ /3). Boundary conditions for the exterior 
of D are p^0as£^+oo, and boundary conditions at £ = £o ar e obtained 
by continuity of p at the boundary. Since for jammed states p = 1 in the 
interior, the continuity condition gives p = 1 on the boundary. Therefore, 
separable solutions must take the form p(£,?7) = F(£/f3)/F(£ /f3). Since the 
system of coordinates is orthogonal and p = 1 in the interior of D from ()34|) . 
we can compute the amplitude of (^-function for density p at the boundary 
£ — £oj which comes out to be —(3F\^q/(3)/F{^q/(3). Thus, the exact solution 
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for the case when separation of variables of Helmholtz equation is possible 
and the boundary is given by dD = {£ = £ } is 




(£, rj) e intD 
(£, rj) G extD 



(35) 



(36) 



' 1, (£, rj) e intD 

F(£/P)/F(Z /P) 7 (£,7i)eextD 

To be specific, in Table 1 we enumerate all cases for which separation of vari- 
ables in Laplace/Helmholtz equation is possible. As it is known, Helmholtz 
equation in two dimensions is separable in four coordinates only: Carte- 
sian (which we shall not consider here), cylindrical, elliptic cylindrical and 
parabolic cylindrical [2B]- Each of these cases selects a particular shape of 
the patch, as well as a particular form of the function F(£). 



Coordinates 


Shape of D 


Functional form of F(£) 


Cylindrical 


Circle 


Bessel function Ko(r) 


Elliptic cylindrical 


Ellipse 


Modified Matthew Function 


Elliptic cylindrical 


Hyperbolae 


Matthew Function 


Parabolic cylindrical 


Parabola 


Parabolic Cylinder Function 



Table 1: A summary of exact jammed solutions and corresponding shapes in 
two dimensions 

Of course, similar results may be obtained in three dimensions, where 
eleven coordinate systems exist for which the Helmholtz differential equation 
is separable. We shall not go into details here as, first of all, the generalization 
is straightforward, and, second, we are only interested here in one- and two- 
dimensional self-assembly. Three-dimensional analogues of self-assembly may 
be important for some models arising in mathematical biology which are 
discussed at the end of Sec|3 



6.2 Irregular shapes: asymptotic result 

For arbitrary shape of domain D, it is no longer possible to find an exact 
expression for densities p and p. However, we may find an asymptotic result 
which shows that 5-function density on the boundary must be a slowly- 
varying function of the boundary coordinates in the case when (3 - the range 
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of H - is small compared to the size of the patch. We consider a patch of an 
arbitrary simple-connected shape Q with a smooth boundary and introduce 
the coordinate t] along the boundary of the patch dQ. Considerations of 
one dimensional case indicate that the patch is stable if it is jammed, i.e. 
p = everywhere inside the patch, which leads to equation p — 1 inside the 
patch. For now, we consider H(x) to be inverse Hemholtz operator, thus p 
satisfies ()34|) . Thus, inside the patch, we necessarily have p — 1. The decay 
of ~p close to the boundary should be compensated by the presence of delta- 
function concentration at the boundary. Let us assume that the boundary is 
smooth. If £ is the coordinate locally perpendicular to the boundary and the 
boundary is at £ = £cb then we assume that the density near the boundary 
has the form: 



where ind^(r) is the indicator function which equals unity, if r e Q and 
vanishes otherwise. From the condition p = 1 and p = 1 in the interior of Q 
we obtain an integral equation for the strength f(s): 



If r is farther than (3 away from the boundary, then only the contribution 
from the second integral is relevant. However, due to the extremely short 
range of H, this integral is equal to unity and equation (jHEJ) gives an identity. 
Thus, we only need to investigate equation (jHEJ) in the case when the distance 
between r and the boundary is of order (3 or less. Suppose this distance is 
d/3 with d > being a constant of order unity or less. Since H{r' — r) decays 
rapidly away from r for distances larger than j3, we can approximate the 
slowly varying function f(s') by its value at s, and the integral ()38j) as 



Since we are only interested in the immediate neighborhood of the point r(?7), 
we will now assume that the boundary is locally straight and vertical, and 
the interior of the patch is to the right of the boundary. We introduce the 
local coordinates x = /?£ and y = j3rj centered at the point r, so the boundary 
is at x = d and the x— axis is pointing towards the boundary. Note the exact 
form of the kernel H 



p(£~£o) = /foW£-£o)+indo(r) 



(37) 




(38) 




(39) 




(40) 
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We split up the second integral in (J39|) and perform the integration exactly 
by going to polar coordinates (and skipping algebraic details) : 



/ H(v - r')dr' ~ — / K (\/x 2 + y 2 ) dxdy = (41) 
Jn 2n J x <d V v / 

Let us now compute the first integral in (JSTJj) . We have: 



Thus, equation becomes 

m ^ e -\r-r>\/P + 1 _l_ e -\r-ryp = 1 (42) 

so the answer is simply 

M = P- (43) 

Thus, for a smooth boundary, jammed states are obtained by the same prin- 
ciple as in one dimension: the densely packed state p = ~p — 1 inside the 
domain is surrounded by a 5-function layer of strength (3. At first glance, 
the asymptotic answer (J4H|) seems not correspond to the exact answer (jHB^) 
for circular, elliptic and parabolic shapes. Consider, however, the case of a 
circular patch as an example. Assume that the boundary of the patch being 
at r = r . According to (J35|) . the strength of 5-function on the boundary is 

K' (r /(3) K x (rp/g 
J " P K (r /(3) P K (r /P) 

The ratio of Bessel functions converges rapidly to 1 for (3 — > 0. In fact, 
this ratio is extremely close to 1 already when (3 < r$. Since (3 is assumed 
small, for a patch of any reasonable size (larger than several units of (3), 
the approximate result (|43|) is valid to high accuracy (exponential in (3). In 
general, we can only postulate that the accuracy of (|4^|) should be 0(/? 2 ), 
since we neglected the local curvature of the surface. Introduction of local 
curvature effects will change the asymptotic result (}4*3"|) . but this correction 
is bound to be small (of the order j3 2 ) when (3 — > 0. 
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6.3 Numerical simulation of states in two dimensions 



To illustrate these exact and asymptotic results, we have performed fully 
nonlinear numerical simulations for a — 1, f3 — 0.1 and D = 0.01 in two 
dimensions. The results of this simulations are shown in Fig. 8. gaussian 
radial distribution. Evolution due to (|OI7j) deforms this shape into a flat- 
top elliptical shape, which is reminiscent of the solutions derived in Sec J6. II 
We conjecture that elliptical shapes are more stable than circular shapes, 
although a detailed analysis of the two-dimensional stability and selection 
mechanism has yet to be completed. 

7 Conclusion, open problems and further ap- 
plications 

A new model was proposed and analyzed for the collective aggregation of 
finite-size particles driven by the force of mutual attraction. Starting from 
smooth initial conditons, the solution for the particle density in this model 
was found to collapse into a set of delta- functions (clumps), and the evolution 
equations for the dynamics of these clumps were computed analytically. The 
energy derived for this model is well defined even when density is supported 
on (5-functions. The mechanism for the formation of these (^-function clumps 
is the nonlinear instability governed by the Ricatti equation (J19)) . which 
causes the magnitude of any density maximum to grow without bound in 
finite time. At first sight, it may seem that the emergence of 5-function 
peaks in the solution might be undesirable and perhaps should be avoided. 
However, these ^-functions may be understood as clumps of matter, and 
the model guarantees that any solution eventually ends up as a set of these 
clumps. Subsequently, further collective motion of these clumps may be 
predicted using a (rather small-dimensional) system of ODEs, rather than 
dealing with the full non-local PDEs. The question of how many clumps arise 
from a given initial condition remains to be considered. One may conjecture 
that clump formation is extensive; so that each clump forms from the material 
within the range of the potential $, determined by G. On a longer time scale, 
the clumps themselves continue to aggregate, as determined by the collective 
dynamics (fTHj) of weak solutions (fT2*j) . This clump dynamics is also a gradient 
flow; so that eventually only one clump remains. 

A comparison with point vortex solutions of Euler's equations for ideal 
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hydrodynamics in two-dimensions may be made here. Prediction of the in- 
compressible motion of an ideal fluid is governed by a set of nonlinear PDE 
in which pressure introduces non-locality. A drastic simplification of motion 
occurs, when all the vorticity is concentrated in delta-functions (point vor- 
tices) j^ni- The motion of point vortices lies on a singular invariant manifold: 
if started with a set of point vortices, the fluid structure will remain a set of 
point vortices. However, a smooth initial condition for vorticity does not split 
into point vortices under the Euler motion. In the present model, though, the 
physical attraction drives any initial distribution of density towards a set of 
delta-functions, so one is guaranteed to obtain effectively finite-dimensional 
behavior in the system after a rather short initial time. Of course, this time 
depends of the precise form of the long range attraction. 

Because two scales are present in the smoothing functions H and G, 
the effects of boundary conditions warrant further study. For example, the 
H — > 5 limit should also allow formation of boundary layers. These boundary 
issues were avoided in the present treatment by using periodic boundary 
conditions. However, it would be natural in some physical situations to 
apply, e.g., Neumann boundary conditions to the particle flux J. The issue of 
boundary conditions will arise and must be addressed on an individual basis 
in specific applications of these equations in physics, chemistry, technology 
and biosciences. 

Our approach involving a variational principle, energy and competition 
of length scales is relevant for many areas of science. In particular, we have 
recently learned that our variational approach is applicable to bio-sciences, 
in particular, to the theory of insect swarming. In a recent work, Topaz et 
al. |23 have discovered the formation of isolated patches of matter (swarms) 
in a variant of (jb17|) with a nonlinear instead of nonlocal diffusion. The 
variational energy (jSJ) for their model was 



This case could be considered as an extension of our model for H(x) = 5(x), 



with Q'(p) > 0, in particular, Q(p) = p 2 /2. Alternatively, a variant of this 
nonlinear diffusion model can be derived if we start with /i(p) = 1 — p and 



E\p\ = J P G*p--p 2 



= 1 and 
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D — 0. We can then rewrite (jUJ) as 

dp 

jf- + div (pV$) = div (ppV$) , 

which, again, is a modification of swarming models considered in [H^ with 
non-local diffusion. We emphasize that, in our opinion, nonlocality in diffu- 
sion is advantageous, since it allows for simple generalized solutions in terms 
of constants and ^-functions. These generalized solutions dominate both 
the dynamics and statics of the problem and greatly simplify the analytical 
treatment. 
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Figure Captions 



Figurel. Scanning Electron Microscope (SEM) figures of self-assembly of parti- 
cles in nano-channels, courtesy of S. R. J. Brueck and D. Xia. Note 
fully dense clumps separated by voids, evident in cases B, C and F. 
Cross-width of the channels as well as distance between the channels is 
100 nm. 

Figure2. Numerical simulation demonstrates the emergence of particle clumps, 
showing formation of density peaks in a simulation of the initial value 
problem for equation (jUJ) using smooth initial conditions for density 
with I = 1, G(x) = H(x) = exp(— \x\), /i(p) = 1. The vertical coordi- 
nate represents ~p = H * p, which remains finite even when the density 
forms (^-functions. 

Figure 3. Position of inflection point for simulation shown on Fig. 2. The dashed 
line in the bottom shows the minimum value allowed by the finite res- 
olution of the mesh. In this case particular case, this minimum value 
is equal 0.1 with the length of the interval being 10. 

Figure 4. Evolution of a Gaussian initial condition for p(x,0) with /x(p) = 1 — p 
and H = exp(— \x\/e)/(2e) where e = a/10. The solution quickly forms 
a plateau of maximal possible density (p m ax = 1)- 

Figure 5. Stationary equilibrium solution for a = 1, (3 = 0.1 and L = 1. Top 
p(x) (with representation of 5-functions at x — ±L as vertical red lines). 
Middle: $ = (G * p) (x) for the same solution. Bottom: ~p = H * p for 
the same solution. 

Figure 6. Stationary jammed states for a — 1, f3 — 0.1 and L — 1. Top p(x) 
(again, 5-functions at x — ±L are represented by vertical red lines). 
Middle: the potential $ = (G * p)(x). Bottom: p = H * p. 

Figure 7. Convergence to analytic equilibrium solution (circles) for different times 
(colored solid lines, see legend). In simulation, D = 0.01, a = 1, 
P = 0.1 

Figure 8. Two-dimensional evolution of a gaussian initial profile for consecutive 
times with D = 0.01, a = 1, /3 = 0.1 
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